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The fundamental dislocation processes of glide, climb, and annihilation are studied on diffusive 
time scales within the framework of a continuum field theory, the Phase Field Crystals (PFC) 
model. Glide and climb are examined for single edge dislocations subjected to shear and compressive 
strain, respectively, in a two dimensional hexagonal lattice. It is shown that the natural features 
of these processes are reproduced without any explicit consideration of elasticity theory or ad hoc 
construction of microscopic Peierls potentials. Particular attention is paid to the Peierls barrier for 
dislocation glide/climb and the ensuing dynamic behavior as functions of strain rate, temperature, 
and dislocation density. It is shown that the dynamics are accurately described by simple viscous 
motion equations for an overdamped point mass, where the dislocation mobility is the only adjustable 
parameter. The critical distance for the annihilation of two edge dislocations as a function of 
separation angle is also presented. 



I. INTRODUCTION 

Plastic flow in periodic systems is typically medi- 
ated by the motion of line defects or dislocations. The 
largest challenge in developing a meaningful theory of 
plasticity is often linking the microscopic behavior of in- 
dividual, discrete dislocations to the macroscopic plas- 
tic behavior of the system. In atomic and molecular 
crystals for example, understanding the effect of dislo- 
cations on mesoscopic and macroscopic material prop- 
erties involves the treatment of length and time scales 
that capture the relevant dynamics of individual disloca- 
tions ('--^10~^^s,~10~^m) through those that describe the 
macroscopic response of the material (~10^s,~10~^m). 
An important approach to the problem of spanning this 
large range of scales has been to measure the dynam- 
ics of individual dislocations and/or small numbers of 
interacting dislocations on the shortest time scales from 
Molecular Dynamics (MD) simulations [1I3,|1,I3. These 
results are then used as input into coarse-grained, meso- 
scopic simulations such as Dislocation Dynamics (DD) 
1^, which can provide information on systems with 
large numbers of dislocations under the action of experi- 
mentally accessible strains and strain rates. 

In this study, dislocation dynamics are examined on 
length scales comparable to those encountered in MD 
simulations, but over diffusive time scales and using ex- 
perimentally accessible strain rates. This approach pro- 
vides a single framework that removes the vibrational 
time scales, while all of the relevant length scales can 
potentially be reached with greater computing power or 
with more advanced numerical techniques 0. In addi- 
tion to atomic crystals, the results presented here may 



be interpreted in terms of other periodic systems such 
as Abrikosov vortex lattices in superconductors , mag- 
netic thin films 0, 0| , block copolymers [Tl| , oil- water 
systems containing surfactants [ij, and colloidal crys- 
tals. 

The PFC model describes periodic systems of a con- 
tinuum field nature and naturally incorporates elastic 
and plastic behavior. The details of the model have been 
presented elsewhere ^sj, and only the necessary equa- 
tions will be given here. The dimensionless free energy 
functional is written as 



F ^ dx 
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where p is an order parameter corresponding here to den- 
sity and 
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r is a phcnomenological constant related to temperature. 
The dynamics of p are described by 



dp 
dt 
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where ^ is a gaussian random noise variable 
mn,h)C{f2,t2)) - DV^Sin ~ f2)<5(Ti - T2) and 
D = uksTq^^* / \^) which has been largely neglected in 
this study, as will be discussed in Section III. 

In Section II, the details of how the PFC model is 
adapted to numerical simulation are outlined, and in Sec- 
tion III the simulation results for glide, climb, and anni- 
hilation are presented and analyzed. Section IV includes 
a summary, comparison with other recent phase field sim- 
ulations of dislocations, and discussion of further devel- 
opments. 



II. SIMULATION METHOD 

A. Discretization, Initial Conditions, and 
Boundary Conditions 

Eq. ||2Jl was solved numerically in two dimensions 
using the 'spherical laplacian' approximation for 0| 
and a forward Euler discretization for the time deriva- 
tive. Periodic boundary conditions were applied in all 
directions for glide simulations and mirror boundary con- 
ditions were used perpendicular to the climb direction in 
climb simulations. To create a system with a single edge 
dislocation, an initial condition consisting of a hexag- 
onal one-mode solution for y) was applied with N 
atoms/row in the lower half and N-l-1 atoms/row in the 
upper half. The hexagonal state is expressed analytically 
as 
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p{x,y) = A 



cos (qx) cos{—=) cos 



where 



Po (4) 
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q is the numerically determined equilibrium wavenumber 
for a hexagonal state at a given value of r, and po is 
the average density. In glide simulations, the hexagonal 
state was bounded at its upper and lower edges by a con- 
stant, or liquid, state of width approximately 4aj,, where 
ay is the equilibrium lattice parameter in the y-direction 
(Fig.^. The same approach was used in the climb simu- 
lations, except that the liquid was placed along the sides. 
Before applying strain, all systems were allowed to equili- 
brate until their free energy no longer changed with time. 

At a given value of r, the value of po for the 
hexagonal portion of the simulation was set to fall on 
the phase boundary between the hexagonal and hexago- 
nal/constant coexistence regions. The value of po for the 
liquid portion of the simulation was set to fall on the 
boundary between the hexagonal/constant coexistence 
region and the constant phase region. This was neces- 
sary to make the interfaces between the hexagonal and 
constant phases stable, with no preference toward crys- 
tallization or melting. A drawback is that this makes 
any comparison of results at different r values indirect, 
since po must vary with r. For this reason the boundary 
conditions were changed when the r dependence of the 
dynamics was of interest. Details are discussed in the 
following section. 



B. Strain Application 

Two methods were used to apply strain to the sys- 
tem. In both, p{x, y) was coupled to an external field 
along the outer two rows of particles bounding the liq- 
uid phase on each side of the system. This external field 
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FIG. 1: Schematic of (111) plane in a FCC crystal corre- 
sponding to the 2D system of interest. 



was set to the one-mode solution given in Eq. and 
for glide (climb) was moved in the positive x-direction 
along the lower (left) rows and in the negative a;-direction 
along the upper (right) rows, both at the same constant 
velocity. The particles in the system are energetically 
motivated to follow the motion of these fields, giving the 
effect of a physically applied strain. 

In the first method, which will be called rigid dis- 
placement, Eq. Q was solved in the presence of the 
external fields, but in addition, the particles between the 
external fields were rigidly displaced along with the mo- 
tion of the fields to ensure a linear strain profile across 
the width of the system. In the second method, this rigid 
displacement was not enforced, allowing the strain profile 
to take whichever form the dynamics of Eq. (PJ dictate. 
This method will be called relaxational displacement. In 
Section III, it will be shown that the dynamic behav- 
ior of the dislocations can be significantly influenced by 
which method is used and that the two methods may be 
viewed as limiting cases of rigid and diffusive response. 
From this viewpoint, rigid displacement describes atomic 
crystals and relaxational displacement applies to 'softer' 
systems such as colloidal crystals, superconducting vor- 
tex lattices, magnetic films, oil-water systems containing 
surfactants, and block copolymers. 



C. Symmetries and Time Scales 

The crystalline symmetry here is equivalent to the 
{111} family of planes in a FCC lattice or the {0001} 
family of planes in a HCP lattice, for example. These 
close packed planes and the subsequent glide directions 
compose the primary slip systems in many common types 
of ductile, metallic crystals. Using the FCC lattice as a 
reference, application of shear in this geometry results in 
glide along a (110) direction within a {111} slip plane, as 
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shown in Fig. ^ The directions in a HCP lattice would 
fall in the (1120) family. Climb in this geometry was 
made to occur along a (112) direction. Shear and com- 
pression were also applied over various other rotations as 
will be discussed briefly in the following. 

System sizes ranging from 676 to 56,952 particles 
were examined, and strain rates ranging from 2xl0^^/t 
to Ix 10~^/i were used, where t is the dimensionless time 
introduced in Eq. Q. These strain rates can be ex- 
pressed in physical units by matching the time scales 
of the model to those of typical metals near their re- 
spective melting temperatures. This is done by equating 
vacancy diffusion constants, Dy, which have been cal- 
culated analytically for this model in [23|, and which 
may range from 10~^-10~^^cm^/s for typical metals [2]| . 
Lattice constants, a, must also be equated to return 
to physical units. Using Cu at 1063°C as a reference 
(Dy ~ 10~^CTO^/s, a ~ 3. 61 A), and matching to the 
model at r = —0.8 {Dy = 1.78a^/i), the range of strain 
rates used converts to .09/s-4500/s. Using these same 
parameters, the dislocation velocities observed are on the 
order of 10~^-10"'*m/s, a range well below the acoustic 
limit and accessible by experiment. The dislocation den- 
sities range from approximately lO^^-lO^^/cm^. 



D. Simulation Output: A Preliminary Example 

Before presenting the analysis of all simulation data, 
the output from a single glide simulation will be pre- 
sented to clarify various definitions and results that will 
be of importance in interpreting the data. The collective 
results from all simulations will be analyzed further in 
the following section. 

Four primary types of output were generated in each 
simulation, from which all properties of interest were ex- 
tracted. The variables are the instantaneous position and 
velocity of the dislocation, and the strain and change in 
free energy of the system, all recorded as functions of time 
as shown in Fig. [3 The gray lines represent theoretical 
results which will be presented in Section III. 

The position was determined by locating all maxima 
of p{x, y) (which can be considered the discrete particle 
locations) and counting the number of nearest neighbors 
for each. Any maxima with more or less than six nearest 
neighbors must be near the dislocation core, and by av- 
eraging the positions of all maxima identified in this way, 
an overall dislocation position was inferred. The velocity 
was then calculated from the slope of the position versus 
time data. 

The average shear strain in the system, 7, was mea- 
sured by again locating the peaks in p(x,y), and not- 
ing that in equilibrium, each particle will have another 
particle located a distance of 2ay away in the positive 
y-direction. If this particle is found to be offset some dis- 
tance, dxi, in the x-direction, then the local shear strain 
is equal to dxi/2ay. The average shear strain in the sys- 
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FIG. 2: One set of simulation data (black lines) and the 
corresponding theoretical results (gray lines) for glide. Pa- 
rameters are r = -0.4, (L^, L;,)=(60,56), and 7 = 2 x 10~^/t 
The inset in the upper left corner of the upper left plot shows 
a magnification of the position vs. time data at early times to 
emphasize the stick-slip nature of the motion at low velocities. 



tern is then given by 



N 



2Nay ^ 

^ Z — 1 



dxi 



(6) 



where iV is the number of particles in the system. The 
fourth variable, the average free energy F ^ was simply 
calculated from Eq. at regular intervals of time. 

The Peierls barrier is a measure of the resistance to 
the onset of motion in a periodic system. In these sim- 
ulations, the barrier is defined as the amount of strain 
that has been applied at the instant that the dislocation 
has precessed a distance of one lattice constant. 7p and 
ep will denote the Peierls barriers for glide and climb, 
respectively. For clarity, Fig. |3| shows 7 as a function 
of time for a few different values of 7. The strain cor- 
responding to this definition of 7p is indicated on each 
curve and can be seen to correspond to the point where 
the measured strain begins to deviate from the applied 
strain. The deviation is due to the strain relieved by the 
motion of the dislocation, as will be discussed in the next 
section. 



III. RESULTS AND ANALYSIS 

A. Equilibrium Dislocation Geometry 

Following equilibration as described in the previous 
section, the dislocations were found to reach one of the 
two stable configurations shown in Fig. 01 Which of the 
two configurations is selected depends sensitively on the 
details of the boundary conditions as well as on the sys- 
tem size. Systems larger in the x-direction tend to favor 
Config. 1, and systems larger in the y-direction tend to 
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FIG. 3: Corrected and uncorrected 7 versus time at r = —0.4 
and {Lx, Ly) = (60, 56) for various values of 7. The strain at 
which 7p is defined has been highhghted on each curve. 



FIG. 5: Equihbrium elastic strain, 7, due to an edge dis- 
location plotted as a function of inverse system size in the 
y-direction. Lx was fixed at 56 particles for the data shown. 
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FIG. 4: Stable dislocation configurations: The greyscale rep- 
resents p{x,y) and the particles around the dislocation core 
have been highlighted for clarity. Left: Config. 1. Right: 
Config. 2. 



favor Config. 2, apparently due to the greater strain re- 
lief available at larger extensions. Systems with approx- 
imately equivalent x and y dimensions that were equili- 
brated with thermal noise oscillated between Configs. 1 
and 2, indicating that the two states are approximately 
equivalent energetically. It will be shown in the following 
section that the initial configuration affects "fp but not 
the velocity of the dislocation. 

The average shear strain, 7, in each system was 
measured and the values recorded following equilibration 
have been plotted in Fig. (Sjas a function of where 
Ly is the number of particles in the y-direction. A simple 
analysis reveals that the total 7 due to an edge disloca- 
tion in this geometry is equal to ^/3b/ {4Ly) , where b is 



the burger's vector of the dislocation. This result agrees 
well with the measured values shown in Fig. [3 indicating 
that the measurement technique is reliable. 



B. Glide: Constant Applied Shear Rate Dynamics 

Simulations were conducted using steady shear over 
a range of applied shear rates (7) , temperatures (r) , and 
system sizes {L^, Ly). The dependence of the Peierls bar- 
rier and the velocity vs. 7 behavior on these variables will 
be discussed in the following subsections. 



1. Peierls Barrier for Glide 

To test for finite size effects, 7p was measured as 
a function of system size, or inverse dislocation density. 
Within estimated errors, no change was observed under 
rigid displacement as the system size was increased from 
676 to 56,952 particles. Under relaxational displacement, 
a slight increase with Ly was noted, and is linked to the 
time required for the strain applied at the edges to diffuse 
inward to the dislocation core. Diffusion is fast compared 
to the inverse shear rates required to apply relaxational 
displacement (rows of particles slip relative to each other 
at all but the lowest values of 7), so the increase of 7p 
with Ly cannot be very large. The nonlinear shear profile 
that is produced may exaggerate this lag between the 
applied strain and the strain near the dislocation, but 
the overall effect was nonetheless found to be relatively 
small. 

Next, the barrier was examined as a function of r, 
which is proportional to the distance in temperature from 
Tc- To do this consistently, the boundary conditions were 
changed to mirror rather than periodic at the top and 
bottom, and the constant phase was entirely removed 
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FIG. 6: Temperature dependence of the Peierls strain barrier 
for glide without thermal fluctuations. Data shown is at po = 
0.25 and {Lx, Ly)={56,56) under rigid displacement. 



from the simulation. This made it possible to vary r at a 
single value of po, isolating the temperature dependence 
in a more realistic manner. Results are shown in Fig. 

The decrease in 7p as the melting point is ap- 
proached is expected since the hexagonal phase becomes 
less pronounced near Tc- That is, A decreases with in- 
creasing r according to Eq. (|3J), and even without ther- 
mal fluctuations a distinct temperature dependence is 
produced. This decrease in A corresponds to an increase 
in the width of the dislocation which, according to the 
Peierls-Nabarro model , lowers the Peierls barrier for 
glide. With thermal fluctuations, these results did not 
change significantly, though at low 7, which is where the 
change would be greatest, it was not possible to include 
fluctuations and maintain reasonable computation times. 
Similar linear decreases in 7p as some effective Tc is ap- 
proached have been found in experiment p^llTj and the- 
ory [TMll9| . along with increases in jp with 7 much like 
those shown in Fig. 

At temperatures closer to the melting point (r ~ 
—0.18), the dislocations began to climb at very low 
strains before any glide had occurred. This is the first 
evidence that climb is the dominant process at high tem- 
peratures, as in real crystals. Further evidence will be 
presented with the climb results. 

The dependence of 7p on 7 was also explicitly mea- 
sured (Fig. O. Both methods of shear application result 
in what appears to be a power law increase in 7p as the 
shear rate is increased, where the relaxational displace- 
ment data is nearly linear and the rigid displacement data 
appears to approach a limit 7p at high 7. 

This increase is explained as follows. Extrapolating 
the data to 7 = indicates that in all cases there is 
some small strain jp < 0.5% at which the dislocation 
will glide, given sufficient time. Call this time from when 
7p is reached to when the first glide event actually occurs 



FIG. 7: Measured Peierls strain barrier for glide, yp, as 
a function of applied strain rate for the cases of rigid and 
relaxational displacement. The fits to the relaxational data 
are from Eq. Q and the fits to the rigid data are power laws 
as indicated in the image. The dotted line shows a linear 7 
dependence for reference. Note that 7p is consistently lower 
for Config. 2 than for Config. 1, even with the differing r's 
working toward the opposite effect. 

Athop- In any given simulation, once 7 > 7p, excess 
strain is being applied during the interval At^op which 
makes the observed 7p appear to be larger than jp. In 
this approximation 

7P = 7pW+7At/,op(r) (7) 

which predicts a linear increase in 7p above 7p. It is 
reasonable to expect though, that Athop will decrease at 
higher strain rates due to the additional strain applied 
during the interval. In a first approximation then 

Athopir, 7) = A<Lp(^) - c^irHd, 7)7 (8) 

where At^Qp(r) is the time required to execute a glide 
event at 7 = and a{r) is a coefficient related to the 
magnitude of the additional driving force applied during 
Athop- 7) is a coefficient related to the efficiency 
of strain transfer from the sample edges where strain is 
applied to the dislocation core, as a function of the type of 
displacement d and 7. r](d,'j) = 1 for rigid displacement 
and drops below 1 under relaxational displacement. 
Substituting Eq. © into Eq. gives 

7P - 7p(0 + i^tl^pir) - airMd, 7)7^ (9) 

which now indicates a dependence on 7 that must be fall 
below the initial expected linear trend. And the devia- 
tion from linear will be greatest under rigid displacement, 
since rj is maximum in this case. The data is in agreement 
with this expectation. Using 7p(— 0.4) = 7p(— 0.8) — 
0.0038, Athop{-OA) = 7277i, Athop{-0.8) = 6643i, and 
q;(— 0.4)?7 = a(— 0.8)ry = 1.48 x 10~^t'^ produces reason- 
able fits to the relaxational displacement data shown in 
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Fig. 13 Note that if At^op is negligible, then clearly this 
effect will not be noticeable, but since the dynamics are 
necessicarily diffusive in these simulations, it is reason- 
able to expect some contribution from this effect. 

Accurate fits to the rigid displacement data were 
more difficult to obtain, most likely due to the transition 
to no 7 dependence for large 7. This can be shown by 
studying the evolution of F under an applied shear. In 
[20I , the change in F for a one-mode approximate hexag- 
onal solution under the action of shear was be found by 
minimizing F when p{x, y) is replaced with p{x -\- jy, y). 
The resulting equation, valid for small 7, is 



AFf 



Shear 
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(10) 



In principle, this represents a rigid displacement of p{x, y) 
at infinitely large 7. In this limit, has no explicit 
dependence on 7; 



7p 



(11) 



Dislocation dynamics in soft structures such as col- 
loidal crystals are reasonably expected to correspond to 
the case of relaxational displacement. These systems typ- 
ically exhibit very little rigidity associated with sound 
modes or phonons, thus their relative softness. Con- 
versely, dynamics in atomic crystals are believed to corre- 
spond to the case of rigid displacement at large 7. Atomic 
crystals exhibit a significant rigidity in response to an 
applied shear, which can be reasonably approximated 
by a linear shear profile as is done for rigid displace- 
ment. A more constructive way to model atomic crystals 
would be to explicitly consider a phonon or wave term 
in the dynamics, as is being done by other authors [2^ . 
If such modes were considered, the collective motion of 
particles in response to an applied force would naturally 
be enhanced, more resembling the case of rigid displace- 
ment. It is argued in this sense that the methods of rigid 
displacement in the large 7 region and relaxational dis- 
placement represent limiting cases of response, and that 
a more rigorous description including effective phonon 
dynamics would fall between these limits. 



2. Atomistic Glide Mechanism 

The nature of the dislocation motion in these simu- 
lations (Fig. ISJ is stick-slip at low velocities with a tran- 
sition to a more continuous character at high velocities. 
This is expected, as the lattice barrier leads to thermally 
activated motion when AFshear approximately equals 
^pGhde ^ while at large values of AFshear the barrier 
becomes secondary and the motion assumes a damped 
character. The shear rate dictates the maximum veloc- 
ity and therefore the extent to which the motion becomes 
continuous. Three regimes of motion were observed, with 
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FIG. 8: Atomistic glide mechanism under constant applied 
shear rate (the particles around the dislocation core have been 
highlighted for clarity). From left to right, p{x, y) is shown at 
t=0, 500, 1000, and 1500, corresponding to 7=1%, 2%, 3%, 
and 4%. The arrows indicate relative strain magnitudes and 
directions. 
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selection depending on the ratio 

7 

Vss = — r 
Pdb 

where pd is the dislocation density dictated by the system 
size. The reason for labeling this quantity v^g will soon 
become apparent. 

For large Vss [f^ 0.016a/i), the dislocation quickly 
reaches the overdamped regime and adjacent layers of 
particles begin to slip relative to each other along the 
x-direction before a steady-state velocity is achieved. 
Slipping usually occurs when the strain exceeds ap- 
proximately 20% in rigid displacement or 10-15% in 
relaxational displacement. At moderate values of Vss 



(O.O67/ 



< 



^ 0.016a/t), the dislocation approaches a 



continuous glide motion and eventually reaches a steady- 
state velocity. This velocity can be calculated by equat- 
ing the Orowan equation 



iPlastic = Pdbv 



(13) 



to the applied shear rate, giving the quantity Vss defined 
in Eq. l(T^ . This is the glide velocity required to plas- 
tically relieve strain at exactly the same rate at which 
it is being applied. Fig. |51 shows v^s versus 7 as mea- 
sured from simulations. The measured values follow a 
linear trend as Eq. H12() predicts, with the slopes in good 
agreement with the theoretical values. This again shows 
that the plastic strain relief due to glide is correctly re- 
produced and that the proper steady-states are achieved. 

At low values of Vgs (^S O.O67P), a more surpris- 
ing type of motion occurs in which the dislocation over- 
comes the Peierls barrier, glides a short distance, and 
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FIG. 9: Measured steady-state glide velocities for two system 
sizes. The upper data points are at r = —0.4 and the lower 
data points are at r = —0.8. 



then comes to a stop. The cycle then repeats itself once 
enough strain is re-accumulated to overcome the barrier 
again. This oscillatory motion will occur whenever the 
velocity assumed just above the Peierls barrier is greater 
than the theoretical v^s for the system. The rate at which 
the dislocation glide relieves strain is temporarily greater 
than the applied strain rate, so the energy falls below the 
Peierls barrier and glide is no longer possible until the 
strain energy again increases sufSciently. 



3. Viscous Dynamics 

Empirically, dislocation glide velocity is described by 
the following equation: 



V = VsiTcft/TsY 



(14) 



where Vs is the shear wave velocity. Toff is the effective 
shear stress on the dislocation, Tg is the material stress 
constant, and m is the stress exponent The stress 
exponent has been found to range from less than 1 to 
over 100 in some cases. For typical pure metals such as 
aluminum or copper, m ~ 1-5. These values may change 
significantly depending on temperature, stress range, and 
local defect densities. For example, in iron, m falls into 
one of three regions (77i<l, m—l, m>l) depending on 
the conditions examined As will be shown in this 
subsection, the dislocation velocity was found to be ap- 
proximately linear in both stress and strain (m ~ 1) for 
all parameter ranges studied. This is not unexpected, as 
higher values of m are often attributed to effects such 
as jogs, impurities, and other defects which modify the 
dynamics from those expected for pure, two-dimensional 
crystals. 

The dynamics of a single gliding dislocation are well 
described by the equation of motion for a point mass in 



rricsvit) = Fq- I3v{t) 



(15) 



where rriofj is an effective dislocation mass, is a con- 
stant proportional to 7, and /3 is a damping constant. 

Equations for v{t), x{t), and j{t) can easily be de- 
rived from this starting point, but first the Orowan equa- 
tion will be used to write rricff , Fq^ and (3 in terms of more 
meaningful parameters. It will be shown that the velocity 
is linear in 7, but assuming this for now, one can write 



vi^) = M^{^{t) - 70) 



(16) 



where is the slope, which can be interpreted as an 
effective mobility for glide. Next note that ^{t) is a func- 
tion of the applied strain and the strain relieved by the 
gliding dislocation; 

^{t)=^t~ pdbx{t). (17) 

Substituting Eq. (|17() into Eq. (|16|l and differentiating 
gives 

v{t) = M^-^ - M^pdbv(t) (18) 
and equating terms in Eqs. (|18|l and H15|) shows that 



and 



Fo 



J_ 

rries 



M-fPdb. 



(19) 



(20) 



This analysis indicates that the damping experienced by 
the dislocation is a result of the strain relief connected to 
the glide process and is not directly linked to the dynam- 
ics of Eq. 0. That is, if the second term on the right 
hand side of Eq. Ijl?!) were removed then both the veloc- 
ity and the strain would be linear functions of time, with- 
out any effective damping. Including this term means 
that the effective damping can be controlled by changing 
Pd, with larger values of pd corresponding to increased 
damping. 

Solving Eq. (|15|l in terms of these new parameters, 
and applying the initial conditions v{0) — and x(0) = 
gives 



V{t) = Vssil 



-Myp^bt 



and 



x{t) 



t + 



1 



-Mypdbt 



7 



M-yPdb^ J M^plb^' 

Substituting Eq. into Eq. I^TJ then gives 



lit) 



1 



-M^Pdbt 



) +70- 



(21) 
(22) 

(23) 



Finally, comparing Eqs. and produces the de- 
sired linear relation assumed in Eq. (|16|l and the similar 
relation Vgs = Myjss, where jss — i/M^pdb. The data 
shown in Figs. Illl andl^verifv that these linear relations 
are observed. 
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FIG. 10; Additional comparisons between simulation data 
(black lines) and viscous motion equations (gray lines) for 
glide where where is the only adjustable parameter. From 
left to right in each plot curves are shown for 7 = 2 x 10^^, 
6xl0"^ 2xl0"^ IxW-^, 5xlO"^ and 2 xlO-'^/t exceptor 
the lower right which shows data for 7 = 2 x 10~*, 1.2 x 10~*, 
8 X 10"^, 4 x lO"'^, and 2 x W'^t. The inset in the lower 
right corner of the lower right plot shows data for lower shear 
rates, 7 = 6 x 10"'', 2 x lO"'', and 5 x 10"^/t where Eq. ^ 
begins to fail. In all plots r — —0.8 and {Lx, Ly) — (56,4:6). 
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FIG. 11: Dislocation glide velocity under rigid displacement 
as a function of the measured average shear strain, 7. A num- 
ber of system sizes, temperatures, and shear rates are shown 
to illustrate the uniformity of the dynamics. The heavy line 
is a representative linear fit. Inset: Dislocation glide veloc- 
ity under relaxational displacement at various shear rates and 
two values of r. The strain values are overestimated due to 
the nonlinear shear profile produced by this type of shearing, 
but the slopes are relatively unchanged. The heavy line is the 
same linear fit as in the larger graph. 



In all of these equations, the only adjustable param- 
eter is M-y, the effective mobility of the dislocation. Using 
values of measured from simulations, Fig.|21shows ex- 
cellent agreement between these analytic results and the 
simulation data for one parameter set, and Fig. 1 101 shows 
similar agreement for various other parameter sets. If it 
is assumed that the free energy obeys the relation to 7 
given in Eq. (|10|l . then Eq. H23|l can be substituted into 
Eq. (dni) to give 



^Fshear — 77 

6 



M^Pdb 



(1 



(24) 



which agrees relatively well with the high shear rate data, 
as shown in Fig. ^] The inset in the lower right of Fig. 
IIUI shows how the agreement begins to fail at lower shear 
rates. This anomaly in the low 7 glide data is not fully 
understood. 

It is worth examining the strain dependence of the 
velocity further. In gradient systems, the velocity of finite 
structures is expected to be proportional to the driving 
force applied Fd, which in this case can be interpreted 
as the derivative of the change in free energy due to the 
application of shear; 



Fn = 



dAF. 



Shear 



1- 



(25) 



Additionally, Eqs. H15|l -I|23 () indicate that velocity is in 
general linear in 7 for this type of overdamped system. 
All simulations resulted in approximately linear veloc- 
ity [v) vs. 7 behavior for dislocation glide, as shown 
in Fig. 111! It is important to correct the overall strain 



shown in Fig. ^2 for that relieved by the glide of the dis- 
location (Eq. H13|l ). especially when using small system 
sizes. 

Both methods of displacement produce nearly the 
same value of under all conditions, though it is more 
difficult to determine the local strain around the disloca- 
tion for the case of relaxational displacement, due to the 
nonlinear shear profile. For rigid displacement, the free 
energy follows the expected form {AFshear ^ 7^) and 
the velocity appears to be linear for 7 less than ^ 10%. 
For relaxational displacement, the anomaly in the free 
energy behavior noted above complicates the results, but 
the velocity remains linear in 7 with values of My sim- 
ilar to those found for rigid displacement. An analytic 
calculation of Al-y would complete this analysis of the 
dynamics, but since the simulation results indicate no 
strong dependencies on any variables, Mj = 0.06ax/t is 
believed to be a reasonable estimate for most cases of 
interest. 

Shear was also applied along directions not lying on 
one of the axes of symmetry with predictable results. As 
the angle 9fi is increased (with 0° denoting alignment 
with a symmetry axis), the Peierls barrier grows but the 
slip direction remains along the nearest symmetry axis. 
Once 9ii becomes large enough, approximately 10 — 30° 
depending on the value of r, the dislocation prefers to 
climb rather than glide, with motion in the general di- 
rection of the applied shear. 

A similar analysis to that presented in Eqs. H15(l - 
(j24|l can be applied to the case of constant strain by 
removing the external force from Eq. (|15|l . The resulting 
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FIG. 12: A small sample climb simulation setup where p(x, y) 
has been plotted. The extra row of particles terminating at 
the core of the dislocation has been highhghted. 

equations were also found to agree well with simulation 
data. It is also worth noting that the velocity vs. 7 
behavior is essentially the same as that shown in Fig. 1111 
when the shear condition is one of constant strain. 



C. Climb: Constant Applied Strain Rate Dynamics 

Climb simulations were conducted using steady com- 
pression over a range of parameter values similar to those 
used for glide simulations. Before presenting the results, 
a caveat on this portion of the study is in order. It 
was found that the results varied systematically (i.e. the 
Peierls barrier decreased) with the grid spacing Ax, ap- 
parently due to the decrease in relevant dimensions with 
compression. A grid spacing small enough to overcome 
this effect could not be reached since the time step must 
be dramatically decreased with Ax. But the nature of 
the results and the essential physics remain the same; 
the data is only shifted by this effect. An example of the 
climb simulation geometry is shown in Fig. 1121 



1. Peierls Barrier for Climb 

The dependence of ep on pd is of the same nature as 
that found for 7p. No change was found under rigid dis- 
placement as Lx and Ly were increased, but an increase 
with Lx was observed under relaxational displacement, 
again in proportion to the diffusion time from the edge 
of the sample to the dislocation core. 

The r dependence of ep is shown in Fig. 1131 for 
various strain rates. Comparison with the glide Peierls 
barrier data in Fig. confirms the same general linear 
behavior, ep is quite large at low r but decreases to- 
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FIG. 13: Temperature dependence of the Peierls strain bar- 
rier for chmb without thermal fluctuations. Data shown is at 
po = 0.25 and (ia;, Lj,) = (52,103) under rigid displacement. 

ward Tc such that there is a crossover close to where 
ep becomes less than 7p. Thus climb is predominant at 
high temperatures, in agreement with the accepted phe- 
nomenology . This was also confirmed in the glide sim- 
ulations where climb was found be preferred near Tc, even 
at very low values of applied shear. Note that the data 
shown in Fig. 1 131 was obtained using modified boundary 
conditions of mirror on all sides with no liquid phase. 

Following _20J, the change in F under compression 
can be calculated by substituting p{x/[l + e), y) into Eq. 
l[T|l and minimizing with respect to A. The result is sim- 
ilar to that for shear; 

^Fcom,. = ^e-^ (26) 

In this limit, as was also the case for glide, ep can be 
written in the form 

/ 2Afg'''""-(po,r)" 

^"^V — — • ^''^ 

The strain rate dependence is also similar to that for 
glide, as shown more clearly in Fig. ^1 The results show 
that 7p ~ -^0-30^ which is similar to the dependence 7p ~ 
^0.37 measured for glide at the same r. The absolute 
values of ep are significantly higher than those for glide 
in this case because of the low value of r that was used. 
The same arguments leading to Eq. lO should apply to 
climb as well, since the behavior seems to be essentially 
the same as that found for glide. 

2. Atomistic Climb Mechanism 

Dislocation climb is a nonconservative process. It 
requires either the diffusion of particles away from the 
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FIG. 14: Measured Peierls strain barrier for climb and the 
subsequent mobilities under rigid displacement a.t r — —1.2 
and (L^,L,^)=(52, 166). 



dislocation core or toward it, unlike glide which involves 
only rearrangements of particles around the core. The 
mechanism of climb is shown in Fig. 1151 where in these 
simulations mass diffuses away from the core since the 
strain is applied through compression. 

Again, similar to what was found for glide, the mo- 
tion has a stick-slip character at low velocities and be- 
comes more continuous at higher velocities. The motion 
proceeds by alternating between configurations 1 and 2 
(Fig. 21). Starting from Config. 2, as shown in the upper 
left image of Fig. ^1 the particle marked with an 'X' dif- 
fuses away, leaving the core in Config. 1 as shown in the 
next image. The two particles marked with arrows then 
merge together, returning the dislocation to Config. 2 as 
shown in the subsequent image. The process repeats as 
long as there is sufficient strain energy to maintain mo- 
tion. For climb in the opposite direction, particles diffuse 
in and split rather then diffuse away and merge, respec- 
tively. This merging and splitting of particles may seem 
unphysical, but in a time-averaged sense these motions 
simply represent diffusion of mass away from or toward 
the dislocation core, which is the fundamental limiting 
process in dislocation climb. 



3. Viscous Dynamics 

The dynamics of a single climbing dislocation are 
well described by the same damped equation of motion 
used to describe glide (Eq. ifT^l) '). Again, the only ad- 
justable parameter is M^, the effective mobility for dislo- 
cation climb. Fig.llGlshows the agreement between these 
analytic results and typical sets of simulation data. 

The velocity versus e behavior shown in Fig. 1161 ap- 
pears to be slightly nonlinear, but this is due to the rela- 
tively short range of motion that could be captured with 
computationally tractable system sizes. An approximate 
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FIG. 15; Atomistic climb mechanism under constant ap- 
plied strain rate. From top left to bottom right, p{x,y) is 
shown at t=300, 600, 800, and 900, corresponding to e=2.4%, 
4.8%, 6.4%, and 7.2%. The particles around the dislocation 
core have been highlighted and the rows near the core have 
been labeled for clarity. The particles marked with an 'X' 
are those which diffuse away between subsequent images, and 
those marked with arrows merge together. 



Mj can nonetheless be extracted, and the results indicate 
first of all that the values of are an order of magni- 
tude higher than those measured for {M^ ~ 0.5). The 
slopes of the v versus e curves are much steeper for climb 
than for glide, but at the same time the velocities remain 
zero to much higher strains due to the larger values of 
ep (except near Tc). Also, is not quite as unchang- 
ing as Mj, in that relatively weak, though measurable 
dependencies on r and e were found. The data indicates 
a slight decrease in Afj with increasing r and an increase 
with e that goes like VI (Fig. IT^ . 

To calculate the dynamics of F, Eq. 1)23(1 can be 
substituted into Eq. to give 



AF, 



Comp. 



M^Pdb 



(28) 



which agrees reasonably well with the data shown in Fig. 
[TBI The difference is mostly due to the low value of r 
used, since the one mode approximation loses accuracy 
away from Tc- No anomaly in F like that found in the 
glide data was observed in the climb simulations. All 
curves of the change in F under compression fall onto 
approximately the same curve when plotted versus e. 

Compression was also applied along directions not 
lying on one of the axes of symmetry at r = — 0.8. As 
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FIG. 16: Comparisons between simulation data (black lines) 
and viscous motion equations (gray lines) for climb. From 
left to right in each plot curves are shown for e — 2 x 10~^, 
1.2 X 10"*, 8 X 10"^ 4 X 10"^ 2 x 10"^ and 6 x 10"V*- In 
all plots r = -1.2 and (L^^, Ly) = (52,166). 




FIG. 17: Measured critical radii for annihilation at r = —0.25 
and po = 0.25. The configuration shown in the center is the 
reference dislocation at (0, 0) from which 6q was measured. 
The inset shows a schematic of the expected behavior as the 
temperature is increased above the crossover r at which climb 
becomes dominant. 



the angle Or is increased, the dislocation first glides some 
distance proportional to Or in a direction along the near- 
est symmetry axis. Then climb begins along the same 
lattice direction as in the unrotated case, with the value 
of ep increasing only slightly with Or. Nearer Tc it would 
be reasonable to expect less tendency toward the initial 
gliding, as climb becomes the preferred type of motion. 
Generally speaking, the application of strain along irreg- 
ular directions relative to the lattice symmetry results in 
a mixed motion of glide and climb. 



D. Annihilation 

Annihilation occurs when two dislocations having 
opposite burger's vectors merge and eliminate each other. 
There exists a critical separation, Tc, at a given angle, Oq, 
below which annihilation will occur, and this separation 
is in principle a function of the crystal symmetry, type 
of dislocation, temperature, relative velocity, and the lo- 
cal strain field. Results were obtained here for the static 
case {v = 0) at a single temperature and under no applied 
strain, for two perfect edge dislocations. 

Consider one dislocation at some location (0, 0) and 
another at {dx,dy) with opposite burger's vector. In ra- 
dial coordinates the separation can be expressed in terms 
of a distance rg and an angle 6*0. At 6*0 = 0°, annihilation 
occurs by pure glide, and as 0^ is increased a mixed mo- 
tion of glide and climb is required, until = 90° where 
annihilation occurs by pure climb, was determined as a 
function of by increasing the initial separation until an- 
nihilation no longer occurred. Periodic boundary condi- 
tions were used in all directions and the parameters cho- 
sen were po = 0.25, r = —0.25, and {Lx,Ly) = (56,43). 
The equilibrium wavenumber at this po and r would re- 
quire 56.5 particles in the x— direction, so placing 56/row 
in the bottom half and 57/row in the top half produces 



a dislocation with no preset bias toward climb in either 
direction. The results are shown in Fig. 1171 

Despite the unbiasing, rc is asymmetric with a pref- 
erence toward climb in the — y direction. This is ap- 
parently a consequence of the asymmetry of the strain 
field across the x— axis of the dislocation core, where 
there is an enhancement of strain in the lower half-plane. 
The particle positions around the core clearly reflect this 
asymmetry. Note that the details of the strain field will 
be slightly different for a dislocation in Config. 1, but 
the same argument should nonetheless hold. 

The elliptical shape oirc{Oo) is expected since 7p < 
Ep for this parameter set. As r is increased, eventually 
Ep < 7p I a-iid the primary axis of the ellipse should coin- 
cide with the y— axis (climb axis), becoming more ellipti- 
cal as Tc is approached. This expected behavior is shown 
schematically in the inset of Fig. El Moving from the 
inner to the outer ellipse corresponds to increasing r. 

Extending the elliptical approximation and assum- 
ing that Tc is directly proportional to the Peierls strain, 
one can write a temperature dependent equation for rc\ 



JA^^ sin^ 6*0 -I- ^2 cos2 Oq 

where and A^ are the slopes of the Peierls strain versus 
r curves for glide and climb respectively. 

IV. CONCLUSION 

Three fundamental dislocation processes have been 
numerically examined in idealized two dimensional set- 
tings using a phenomenological PFC model. The dif- 
fusive dynamics were measured over a range of tempera- 
tures, dislocation densities, and experimentally accessible 
strain rates. In equilibrium, two stable edge dislocation 
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configurations were found to exist, with one resulting in 
a slightly lower Peierls barrier for glide than the other. 
The Peierls barriers for glide and climb, 7p and ep re- 
spectively, were found to have little or no dependence 
on dislocation density, and both showed approximately 
linear decreases with increasing temperature (in the ab- 
sence of thermal fluctuations). Near Tc, ep < 7p verify- 
ing the expectation that climb is dominant at high tem- 
peratures. A crossover temperature was identified below 
which 7p < ep and glide becomes the preferred type 
of motion. Both strain barriers also showed essentially 
power law increases with the applied strain rate, where 
the exponents are similar for glide and climb at equal 
r's. Under relaxational displacement (no phonons), 7p 
is nearly linear in 7 and goes as Eq. (jnj, while under 
rigid displacement at high strain rates (strong phonons) 
the deviation from linear is much greater (7P ~ 
with relatively little change in the barrier strain at high 
strain rates. Physical arguments and some mathematical 
arguments were given for all of these behaviors. 

Rigid displacement with 7, e > 5 x 10^^/t most ac- 
curately reproduces the rigid behavior of a real crystal. 
A more rigorous PFC model, derived from microscopies, 
that will include a wave term to simulate phonon dynam- 
ics is currently being developed . Based on the results 
presented here, it is expected that this model will pro- 
duce dynamics that fall between the limits of relaxational 
and rigid response. At strain rates below approximately 
2 X 10~''/i, the two methods of displacement are essen- 
tially equivalent. 

The motion of a gliding or climbing edge dislocation 
was found to be stick-slip in character at low velocities 
and nearly continuous at high velocities. Three possi- 
ble regimes of motion were observed for glide, depending 
on the expected steady-state velocity of the dislocation 
defined in Eq. H12() . These involve an oscillatory glide, a 
steady-state glide, and slipping rows of particles, in order 
of increasing Vss ■ 

A simple viscous dynamic model has been formu- 
lated to describe the results obtained for gliding and 
climbing dislocations, where the only adjustable parame- 
ter is My or Me. Excellent agreement is obtained between 
these equations and the simulation results, both of which 
indicate that velocity is linear in strain for both glide and 



climb. The slope of the v versus 7 curve for glide, M^ 
was found to be nearly unchanging across all parameter 
ranges. The slope for climb, M^, which is an order of 
magnitude greater than M^, was found to increase ap- 
proximately as 

A critical distance for the annihilation of two edge 
dislocations was also measured, and an asymmetry with 
preference toward annihilation in the —y direction was 
found. Tc^Oq) approximately takes the form of an ellipse 
whose major axis is predicted to be along the glide direc- 
tion at low temperatures and along the climb direction 
at high temperatures. 

Other phase-field models have recently been used to 
study dislocation dynamics 18, 23, 24,, ^J. These ap- 
proaches differ from the PFC method in that they do 
not naturally contain atomistic detail. The domains in 
these models typically differentiate dislocation loops and 
the interfaces represent dislocation lines. Coarsening of 
large arrays of lines etc. can be efficiently studied, but 
atomistic detail is either lost or must be explicitly added 
through postulated Peierls potentials. The relevant equa- 
tions of elasticity must also be rigorously applied, unlike 
in the PFC model which naturally exhibits elastic behav- 
ior as well as Peierls potentials. 

Other phenomena relevant to dislocation dynamics, 
such as obstacle and impurity effects, could be studied 
with a similar approach, and more complicated dynam- 
ics involving screw dislocations, dislocation loops, mul- 
tiplication processes, etc. could be examined in three 
dimensional simulations. Alternatively, the two dimen- 
sional model could provide interesting insights into the 
problem of dislocation-mediated melting in two dimen- 
sions. 
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